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USE  OF  THE  TENSOR  PRODUCT  FOR  NUMERICAL  WEATHER 
PREDICTION  BY  THE  FINITE  ELEMENT  METHOD  -  PART  1. 


Introduction 

In  Ref.  1  Hinsman  has  developed  a  Finite  Element  (FE) 
program  for  Numerical  Weather  Prediction  applications.  The 
grid  employed  is  rectangular  with  nodes  at  the  intersections 
of  north- south  and  east -west  lines.  It  was  shown  by  Stani- 
forth  and  Mitchell  (Ref.  2)  that  the  coefficient  matrices 
for  such  a  grid  could  be  expressed  as  tensor  products.  In 
these  products  the  factors  are  matrices  which  depend  solely 
on  grid  spacing  in  the  two  orthogonal  directions.  This 
report  deals  with  the  coefficient  matrix  called  the  "mass" 
matrix  in  FE  parlance.  (In  Refs.  3  and  4  applications  of 
the  tensor  product  resolution  to  the  FE  "stiffness"  matrix 
are  considered.)  The  theory  which  underlies  the  economical 
computational  scheme  based  on  the  mass  matrix  resolution  is 
first  presented.  Next,  the  number  of  floating  point  opera¬ 
tions  and  the  number  of  storage  locations  needed  for  the 
coefficient  matrix  of  this  scheme  are  compared  with  those 
required  by  other  better-known  algorithms.  A  set  of  FORTRAN 
subroutines  for  implementing  the  tensor  product  scheme 
(TENSOR)  is  given  in  Appendix  B. 


Theory 


Consider  the  grid  shown  in  Fig. 


n  columns 


1. 


There  are  n  spaces 


(Note  the  cyclic 
boundary  condition 
in  the  E-W  direction.) 


Fig.  1.  Node  numbering  and  spacing. 


along  each  of  e  grid  lines  in  the  east-west  direction.  Node 
numbering  is  from  west  to  east  along  successive  grid  lines, 
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beginning  in  the  southwest  comer.  There  is  a  cyclic  bound¬ 
ary  condition  in  the  east -west  direction  so  that  the  node 
number  appearing  at  the  beginning  of  each  horizontal  row  is 
repeated  at  the  end.  Spacings  of  the  horizontal  and  the 
vertical  grid  lines  are  not  necessarily  uniform. 

The  computational  problem  addressed  here  is  the  solution 
of  the  equation 

Mw  =  v  <1> 

where  M  (t  e  "mass"  matrix)  is  a  square,  symmetric  matrix  of 
size  ne  and  w  and  v  are  column  vectors  of  height  ne.  M  and 
v  are  input  quantities  and  w  is  sought.  The  tensor  product 
representation  of  M  is 

M  =  MB  *  MA  <2> 

where  MB  is  a  square,  symmetric,  tridiagonal  matrix  of  size 
e  and  MA  is  also  square,  symmetric  and  of  size  n.  MA  is 
tridiagonal  except  for  nonzero  elements  in  upper  right  and 
lower  left  comers.  MB  depends  solely  on  the  north- south 
node  spacing  b  and  MA  depends  upon  the  east -west  node  spac¬ 
ing  a  .  The  asterisk  (*)  denotes  the  tensor  product. 
Explicit  expressions  for  matrices  MA,  MB,  and  the  tensor 
product  are  given  in  Appendix  A. 

Let  MB  be  represented  as  (e  =  3) 

mbn  mb  12  0 

MB  *  mb  21  mb  2  2  mb  2  3  <3> 

0  mb32  mbs j 

If  we  partition  w  and  v  into  e  n  x  1  subvectors  so  that 


wi  fa 

I11  V  *  vV" 

"inj  L  ir 


we  may  use  <2>,  <3>  and  <4>  to  rewrite  <1>  as 

mbit  MA  Wj  +  mbn  MA  Wjj  *  Vj 

mbji  MA  Wj  +  mbia  MA  Wjj  +  rabu  MA  ■  Vjj 
mbit  MA  Wjj  +  mbss  MA  Wjjj  »  Vjjj 
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Define 


W  -  <Wj  wn  wIH> 


<6> 


and 

It  is  easy 
to 


V  .  <v,  VU  v,„> 
to  verify  that  the  equations 

MA  W  MB  =  V 


<5>  are  equivalent 
<7> 


Solution  of  <7>  can  be  accomplished  by  standard  Gaussian 


elimination  procedures.  Specifically,  the  following  steps 
are  required. 


(1)  LDLT  factoring  of 

(2)  Forward  reduction 


ii) 


hand  side  vectors 
LDLT  factoring  of 
Forward  reduction 
hand  side  vectors. 


MA  (n  x  n) . 

and  back-substitution 

MB  (e  x  e ) . 

and  back-substitution 


for  e  right - 
for  n  right - 


This  entire  process  is  economical  of  both  storage  and  arith¬ 
metic  operations  because  of  the  tridiagonal  structure  of  MA 
and  MB . 


Boundary  Conditions 

The  cyclic  boundary  condition  is  implemented  by  repeating 
the  node  numbers  of  the  western  boundary  on  the  eastern 
boundary  as  shown  in  Fig.  1.  As  already  noted,  this 
accounts  for  nonzero  entries  in  the  upper  right-hand  and 
lower  left-hand  corners  of  MA. 

It  is  sometimes  required  to  impose  a  Dirichlet  boundary 
condition  on  the  southern  and  northern  boundaries  of  the 
region.  Specifically,  the  subvectors  w  and  w  of  the  solu¬ 
tion  vector  w  are  prescribed.  To  implement  this  boundary 
condition  the  following  modifications  to  the  standard 
solution  procedure  are  required. 

In  the  n  x  e  matrix  V  on  the  right-hand  side  of  <7>  the 
first  and  last  columns  are  replaced  by  the  prescribed  bound¬ 
ary  values  of  w,  i.e.,  put  Vj  =  Wj  and  vg  =  wg.  Let 
X  =  W  MB  and  solve  the  system 

MA  X  =  V  <8> 

processing  successive  columns  of  V  in  standard  fashion,  but 
omitting  the  first  and  last  columns.  The  reduced  problem 


now  takes  the  form  W  MB  =  X  and  the  first  and  last  columns 
of  X  are  Wj  and  we ,  respectively.  Transposing  both  sides  of 
this  equation  gives 

MB  WT  =  XT  <9> 

where  WT  and  XT  are  the  respective  transposes  of  W  and  X 
(recall  that  MB  is  symmetric  and  is  thus  not  altered  on 
transposition).  Since  the  first  and  last  rows  of  WT  are 
known,  the  corresponding  scalar  equations  are  not  needed. 
Accordingly,  we  form  MB1  by  deleting  the  first  and  last  rows 
of  MB.  We  also  reduce  XT  to  XT1  by  omitting  the  first  and 
last  rows.  This  leaves  the  result 

MB1  WT  =  XT1  <10> 

or,  in  extenso,  this  takes  the  form  (for  n  =  3,  e  =  5) 

m 

mb 21  mb 2 2  mb  2  a  0  0 

0  mb  3 2  mb  3 3  mb  3 4  0 

0  0  mb  %  3  mb  4 4  mb  4 5 

(In  WT  and  XT1  the  elements  denoted  by  "k"  are  known  and 
those  denoted  by  "u”  are  unknown. )  This  equation  may  be  put 
in  standard  form  by  first  altering  the  first  row  of  XT1  by 
subtracting  mb2i  times  the  corresponding  entries  in  the 
first  row  of  WT  and  altering  the  last  row  of  XT1  by  sub¬ 
tracting  mbi,s  times  the  corresponding  entries  in  the  last 
row  of  WT.  Calling  the  new  right  hand  side  XT2  and  forming 
MB2  from  MB1  by  discarding  the  first  and  last  columns  and 
forming  WT1  by  discarding  the  first  and  last  rows  of  WT,  the 
result  is 

MB2  WT1  =  XT 2  <11> 

Solution  of  <11>  is  carried  out  by  LDLT  factoring  of  MB2, 
followed  by  forward  reduction  and  back  substitution. 


Floating  Point  Operations  and  Matrix  Storage  Requirements 


Presented  here  is  a  comparison  of  floating  point  opera¬ 
tion  counts  and  matrix  storage  requirements  for  the  tensor 
product  scheme  and  three  widely-used  solution  algorithms  for 

6 


solving  <7>  or  its  equivalent  <1>.  Three  of  the  schemes 
take  advantage  of  symmetry  of  the  coefficient  matrices  and 
store  only  elements  on  or  above  the  principal  diagonal.  One 
of  these,  the  "band  solver"  (BAND),  places  these  elements  in 
a  rectangular  matrix  ne  x  r,  where  r  is  the  maximum  row 
length  of  the  upper  triangle  of  M.  The  "sky-line  solver" 
(SKY)  further  economizes  by  storing  only  that  part  of  the 
upper  triangle  beginning  at  the  diagonal  and  extending  to 
the  topmost  nonzero  element  of  each  column.  These  subvec¬ 
tors  are  assembled  into  a  single  vector.  This  scheme 
requires  an  additional  integer  address  vector  of  length 
ne  ♦  1.  The  remaining  algorithm,  "successive  over- relax¬ 
ation"  (SOR),  is  iterative  rather  than  direct. 

In  most  applications  of  the  direct  solvers  the  number  of 
floating  point  operations  required  to  factor  the  coefficient 
matrix  into  LDLT  form  is  much  greater  than  those  required  to 
complete  the  process  of  finding  a  single  solution  vector  w 
corresponding  to  a  given  right-hand  side  vector  v  (forward 
reduction  and  back- substitution) .  In  the  present  applica¬ 
tion,  however,  the  latter  solution  process  must  be  carried 
out  17  times  for  each  time  step,  so  that  the  LDLT  factoring 
makes  a  negligible  contribution  to  the  total  computational 
expenditure.  Accordingly,  the  operations  required  for  fac¬ 
toring  are  not  included  in  the  tabulation  below. 

In  the  following  table  the  results  given  for  the  number 
of  floating  point  operations  are  given  in  terms  of  the  grid 
parameters  n  and  e  (defined  in  Fig.  1).  One  multiplication 
(or  one  division)  plus  one  addition  (or  one  subtraction)  is 
counted  as  one  operation.  Exact  results  for  these  operation 
counts  would  take  the  form  of  a  polynomial  in  n  and  e.  Only 
the  highest  degree  terms  are  given  in  the  table.  Since  it 
is  not  possible  to  predict  the  number  of  iterations  per 
solution  when  using  SOR,  the  operation  count  given  for  that 
algorithm  is  for  a  single  iteration.  Also,  since  the  number 
of  storage  locations  required  for  SOR  coefficient  matrices 
is  highly  grid - dependent ,  no  such  entry  is  given  for  SOR. 


1 

Algorithm 

Number  of  Operations 
per  Solution 

Number  of  Storage  Locations 

for  Coefficient  Matrices 

SOR 

10  en  (1) 

(2) 

SKY 

2  en2 

en2 

BAND 

4  en2 

2  en2 

TENSOR 

8  en 

3  n  +  4  e 

Notes:  1.  Number  of  operations  per  iteration. 

2.  Number  of  storage  locations  is  grid-dependent. 


Conclusion 

Close  comparison  of  operation  counts  and  storage  require¬ 
ments  leads  to  the  conclusion  that  the  TENSOR  algorithm  is 
clearly  superior  to  the  SKY  and  BAND  algorithms.  The  com¬ 
parison  with  SOR  is  not  as  clear-cut.  Considering,  however, 
that  the  operation  count  for  SOR  is  for  only  one  iteration, 
there  really  seems  to  be  little  doubt  that  TENSOR  is  the 
method  of  choice. 
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APPENDIX  A 

MATRICES  MA,  MB,  AND  THE  TENSOR  PRODUCT 


Symbols  a..  and  b.  which  appear  in  MA  and  MB  are  defined 
in  Fig.  1. 


2(34+81)  3i  0  84 

^  _  1  ai  2 ( a i +  a2)  a2  0 

^  0  82  2(82+33)  83 

(n  *4)  a4  0  83  2(a  3+  a4 ) 


HB  s  |  bx  2(bx+  b2 )  b2 

(..3)  L°  b>  2b: 


The  tensor  product  of  matrices  C  and  D  may  be  represented 
in  block  partition  form  as 

C11 0  Cj 2  D  Ci 3  D 
C*D  *  C21 0  C22  D  c2 3  D 

.c 3 1  D  C32  D  C33  D, 


where  the  c. •  are  the  elements  of  C.  Note  that,  if  C  and  D 

*  J 

have  dimensions  r  x  s  and  t  x  u,  respectively,  the  tensor 
product  has  dimensions  rt  x  su. 


10 


APPENDIX  B 

FORTRAN  PROGRAM  LISTINGS 


FORTRAN  programs  for  the  implementation  of  TENSOR  are 
listed  here.  They  appear  in  the  form  of  subroutines  AMTRX2, 
FACTOR,  BACKA,  and  BACKB  within  the  test  program  GAUSS3. 
(The  subroutines  FACTOR,  BACKA,  and  BACKB  are  adapted  from 
subroutine  COLSOL  of  Ref.  5.)  Also  included  are  G0G3,  an 
Exec  used  to  execute  GAUSS2  -  a  dimensioned  version  of 
GAUSS3,  and  CDIM,  an  Xedit  program  used  to  enter  the  dimen¬ 
sions  . 


V 

V 


Listing:  GAUSS3  FORTRAN 

C  MAIN  PROGRAM  MASS  MATRIX  USING  TENSOR  PRODUCT  FACTORS 
C 

C  THIS  PROGRAM  IS  DESIGNED  TO  TEST  THE  SCHEME  (TENSOR) 

C  WHICH  RESOLVES  THE  MASS  MATRIX  INTO  A  TENSOR  PRODUCT  IN 

C  ORDER  TO  SOLVE  THE  SYSTEM  OF  EQUATIONS  M  w  =  v.  THE 

C  SUBROUTINES  MAY  BE  INSERTED  IN  THE  PROGRAM  DEVISED  BY 
C  HINSMAN. 

C 

IMPLICIT  REAL*8(A-H.O-Z) 

COMMON/ CM1A/ NLAT, NLONG 
COMMON/CM8 /AiZl) .bT? 12  ,  % 

COMMON  AG (ZB ) ,BG(ZC) ,GA(ZK) ,GB1(ZL) ,GB2 (ZL) ,MA(ZM) , 


1000 


1MB (ZN 
DIMEN 
READ? 
LATX= 
WRITE 
FORMA 


NSION  V(ZP) 

(5,  *)NLONG ,  NLAT 
=NLAT+1  v 

AT^/?^°  MASS  MA 


MASS  MATRIX  -  TENSOR  PRODUCT  RESOLUTION’ 


503 

500 

1001 


FORMAT (7 
FORMAT r 


A  * 

NLONG 


, (24F3 .  0  )  ) 
il24F3t0jj 


NLAT  = ’  ,13.  /) 


ii  *  m  X  \  MUViW  ■  X  W  4  Ufm  x  —  t  •  /  I 

CONSTRUCT  FACTORS,  GA,  GB1 ,  AND  GB2,  OF  MASS  MATRIX 

/I  1  1  T  1  «m«k  99  ’  v 


CALL  AMTRX2 
WRITE?6,501)AG 

rntJMAT  /  7  »  '  a, 


1002 

1004 

1005 


FORMAT (/ 
WRITE? 6, 
FORMAT?; 
WRITE? 6, 
2  FORMAi;?; 
,  WRITE? 6. 

4  FORMAT?; 
WRITE?6. 

5  FORMA] ;  (  ; 
WRITE? 6, 
WRITE  6’ 
CU=GBI?3 
CL=GB1?2 
K= (NLAT- 

IF  NDIR>0. 
NORTHERN  An 


l  1902,1?;  •<12F41)> 

P,„  'gA  ./.(3X.6F7.3)) 
f  jl(J95 JGB?1'  ’  /  '(3X'6f7-3)) 

^  (6  3  «iBZ  •/'(3x-6F7-3)) 

6!l006)MB 


6  1094)< 

¥)< 

?2^LATX 

t-iT*nl( 

,  THERE 


ATX- 1) 

*NLONG 

ERE  IS  A  DIRICHLET  BOUNDARY  CONDITION  ON 
SOUTHERN  BOUNDARIES. 


READ (5 , * j  NDIR.V 
WRITE?  4 ,5Q2jNf>IR,  V. 

502  FORMAT? ;,T  NDIR  =  * ,11,/.’ 
C  PERFORM  LDLf  FACTORING  OF  gA, 
CALL  FACTOrTgA.MA, NLONG) 
CALL  FACTOR  GB1, MB, LATX) 
CALL  FACTOR CGB2, MB, LATX) 
WRITE (6, 1002 )GA 


V:  ' ,6F8.2,/ , (4X.6F8.2)) 
GB1,  AND  G62 


WRITE (6, 1002 )GA 
WRITE ? 6 ’ 1004 )GB1 
WRITE? 6 t 1005 )GB2 

PERFORM  FORWARD  REDUCTION  AND  BACK- SUBSTITUTION  USING 
FACTORS  OF  GA, 

CALL  BACKA (GA ,  V , MA . NDIR ) 

DIRICHLET  BOUNDARY  CONDITION  ON  NORTH  AND  SOUTH 
BOUNDARIES  ? 

if7ndir.gt.q)go  TO  3 
WRITE?6 , 5 10 JV 

PERFORM  FORWARD  REDUCTION  AND  BACK- SUBSTITUTION  USING 
FACTORS  OF  GB1 

CALLqB£CKB ( GB 1 , V , MB , NDIR ) 

CORRECT  RIGHT-HAND  SIDE  FOR  DIRICHLET  CONDITION 
DO  2  J=l, NLONG 

V?J>NLONd)=V(J^NLpNG)-CU*Y(J)  % 

V? J +K ) = V? J >KJ - CL*V ( J ♦ NLAT*NLONG ) 


DO  2  J=l, NLONG 

V(  J*NLON<5 )  =  V  ( Ji- NLONG  )  -  CU* 

SiiiIldio^-CL^"J*NLAT 

PERFORM  FORWARD  REDUCTION  AN 
FACTORS  OF  GB2 

CALL  BACKB ( GB2 , V , MB , NDIR ) 


REDUCTION  AND  BACK- SUBSTITUTION  USING 


A  SECOND  SOLUTION 


6  WRITE (6, 5 10 )V  . 

510  FORMAT?/*  V:  ’ . 6F8 . 2 , / , (4X, 6F8 . 2) ) 

C  READ  A  NEW RlGHT-HAND  SIDE  AND  PERFORM 
READ (5  *)NDIR,V 
write?  6  • 502  TnRir  ,  V 
CALL  BACKA(GA,V,MA,NDIR) 
ifTndir.gt.qjgo  TO  7 
WRITE? 6, 5 10)V 
CALL  BA<iKB?GBl,V,MB,NDIR) 

GO  TO  16 

7  DO  5  J=1 .NLONG 

VT J+NL0N6 ) = V ( J+NLONG ) - CU*V (J ) 

5  V? J+K)=V?J+Kj-CL*V( J+NLAT*NLONG) 

WRITE!  6, 5 10W 
CALL  BA<iKB?GB2,V,MB,NDIR) 

WRITE?6 , 5 10 ) V 
FORMAT?/,;  MA:  , 2X, 3613 ) 

FORMAT (  /  ,  MB:  ,2X,36I3) 

STOP 

C  <r***5S5**'*'***'***Tfc****'****'**'*****'**'*'5V"**"<f ********************* 
SUBROUTINE  FACTOR ( A, MAXA.NN) 

C  . 

C  . 

C  . 

C  . 

C  . 

C  .  NN 
C  .  NWK 

C  .  -  -  OUTPUT 

C  .  A (NWK) 

C  • 
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1003 

1006 


NPUT  VARIABLES  -  - 

=  STIFFNESS  MATRIX  STORED  IN  COMPACTED  FORM 
NP)  =  VECTOR  CONTAINING  ADDRESSES  OF  DIAGONAL 
ELEMENTS  OF  STIFFNESS  MATRIX  IN  A 
NUMBER  OF  EQUATIONS 
NUMBER  OF  ELEMENTS  BELOW  SKYLINE 

=  D  AND  L  -  FACTORS  OF  STIFFNESS  MATRIX 


C 

C 

& 


50 


60 


70 

80 

90 


100 


looo 


IMPLICIT  REAL*8 (A-H.O-Z) 

DIMENSION  A(1),MAXaU) 

PERFORM  L*D*LT  FACTORIZATION  OF  STIFFNESS  MATRIX 

DO  140  N= 1 , NN 
KN=MAXA(N) 

KL=KN+1 

KU=MAXA(N+1)-1 
KH=KU-KL 
IF(KH) 110 ,90,50 

IC=0 

KLT=KU 

DO  80  J=1,KH 

IC=IC+1 

KLT=KLT-1 

KI=MAXA(K)  % 

ND=MAXA?K+ 1 ) -KI - 1 
IF (ND) 80, 80 ,60 
KK=MINO(IC,ND) 

Cs0 

DO  70  L=1.KK  , 

C=C+A(KI+Rj*A(KLT+L) 

A(KLT)=A(KLT)-C 

K=K+1 

K=N 

B=0 . 

DO  100  KK=KL,KU 
K=K- 1  x 
KI=MAXA(K)  % 

A(KK)=C,  ' 

A!KN)=A(KN)-B 

IF?AiKN JJ 120.120,140  v 

WRITE(IOUT,2600)R,A(KN) 

STOP 

CONTINUE  . 

FORMAT?//,'  STOP  -  STIFFNESS  MATRIX  NOT  POSITIVE 
pEFJNITE,  ,^^TNONfOSITiyEx PIVOT  FOR  EQUATION 

rRtxw  ’  0T 

END 


.E20.12) 
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q  a********************************************************* 

0  ****22§§22XX35*  S^iiXi54?4J?il2  ********************** 

C  THIS  SUBROUTINE  PERFORMS  THE  FORWARD  REDUCTION  AND  BACK- 
C  SUBSTITUTION  USING  THE  FACTORS  OF  GA 
C 

IMPLICIT  REAL*8(A-H.O-Z) 

COMMON/CM 1A/NLAT  jNLONG 
DIMENSION  A(I)VVU)  }MAXA(1) 

C  REDUCE  RIGHT-HAND- SIDE  LOAD  VECTOR 


JMIN=L 
LATX=NLAT+ 1 
JMAX=LATX 

IS  THERE  A  DIRICHLET  : 

if7ndir.lt. 1) GO  TO 

SKIP  NORTH  AND  SOUTH 
JMIN=2 
JMAX=NLAT 

)  DO  240  J= JMIN , JMAX 

)  DO  180  n=i,nl6ng 

KL=MAXA (N ) + 1 


HLET  BOUNDARY  CONDITION? 

GO  TO  140 

OUTH  BOUNDARIES 


0  N=1,NL0NG 
XA(N)+1 
XA(N+lj- 1 
-KL) 180, 160, 


KU=MAXA(N+1) - 1 
IF (KU-KL ; 180 , 16 

C=0. 

DO  170  KK=KL,KU 
K=K- 1 


C=C+A(KK)*V(K-*‘ ( J-  l)*NLONG)  . 

V f N+Tj - 1 ) *NLONG )  =  V  (  N+  (J  - 1 ) *NLONG ) - C 
CONTINUE 


BACK- SUBSTITUTE 

DO  200  N=l,NLONG 
K=MAXA(N). 


V (N-^  J - 1 ) *NLONG ) = V ( N+ ( J - 1 ) *NL0NG ) / A ( K ) 

L=2 ,NLONG 
KL=MAXA(N)+ 1 


KU=MAXA(N+ll-l 
IF(KU-KL)230,21 

DO  220  KK=KL,KU 


210,210 


r=k-  1 

220  iV£Kj|j-l)*NLONG)=V(K^(J-l)*NLONG)-A(KK)*V(N^(J-l)* 

230  N=N- 1 ' 

240  CONTINUE 
RETURN 
END 
C 

0  ********************************************************* 


0  ********************************************************* 
C  THIS  SUBROUTINE  PERFORMS  THE  FORWARD  REDUCTION  AND  BACK- 
C  SUBSTITUTION  USING  THE  FACTORS  OF  GB1  OR  GB2 


IMPLICIT 
COMMON/ Cl 
DIMENSIO: 


IT  REAL*8(A-H.O-Z) 
/CM1A/NLAT,NL0NG  ' 

ION  A(l)?V(l),MAXA(l) 


REDUCE  RIGHT-HAND- SIDE  LOAD  VECTOR 


LATX=NLAT+1 
NMIN= 1 
NMAX=LATX 
IS  THERE  A  DIRI 

if7ndir.lt. 1 

SKIP  NORTH  AND 
NMIN=2 


NMIN=2 
NMAX=NLAT 
DO  240  J= 
DO  180  N= 
KL=MAXA(N 


HLET  BOUNDARY  CONDITION? 

GO  TO  50 

OUTH  BOUNDARIES 


J = 1 , NLONG 
ljj=^MIN,NMAX 


IF  l^KU-KL)  180 , 160 , 160 
C=0. 

DO  170  KK=KL,KU 
K=K-1 

C=C+A(KK)*V(J+(K-l)*NLONG) 

V ( J*  ?N - 1 ) *NLONG ) =V(J+(N- 1) *NLONG ) - C 
CONTINUE 

BACK- SUBSTITUTE 

DO  200  N=NMIN,NMAX 

V? JV(N- 1 ) *NLONG ) = V ( J+ ( N- 1 ) *NLONG ) / A (K ) 
LMAX=LATX 

IfTnDIR.lt. 1) GO  TO  205 

LMIN=3 

LMAX=NLAT 

N=LMAX 

DO  230  L=LMIN,LMAX 
kl=maxa7n ) + l 

KU=MAXAjNU)-l 
IF(KU-KL)230 ,210 ,210 

DO  220  KK=KL,KU 


K=K-1 

*0  V(JMK-l)*NLONG)=V(J*(K-l)*NLONG)-A(KK)*V(J+(N-l)* 

IQ  N^N-r 
>0  CONTINUE 
RETURN 
END 

****■*****************************■**■********************•&•* 

****f^5525J5i*I***?*5J************************************ 

THIS  SUBROUTINE  FORMS  THE  MASS  MATRIX  IN  THE  FORM  OF  A 
TENSOR  PRODUCT  OF  THE  GB  MATRIX  AND  THE  GA  MATRIX. 

THE  FIRST  OF  THESE  IS  NLAT  +  1  BY  NLAT  ♦  1. SYMMETRIC, 

AND  TRIDIAGONAL.  THE  SECOND  IS  NLONG  BY  NlLONG, 
SYMMETRIC,  AND  TRIDIAGONAL  EXCEPT  FOR  SINGLE  ELEMENTS 
IN  UPPER  AlGHT  HAND  AND  LOWER  LEFT  HAND  CORNERS.  GB  IS 
STORED  IN  SKYLINE  VECTOR  FORM  (UPPER  TRIANGLE  WITH  SPACE 
FOR  FILL-IN)  AS  GB1  AND  GB2 .  THE  LATTER  VERSION  IMPOSES 
A  DIRICHLET  BOUNDARY  CONDITION  ON  THE  NORTH  AND  SOUTH 
BOUNDARIES.  GA  IS  ALSO  STORED  IN  SKYLINE  VECTOR  FORM. 
INTEGER  ADDRESS  VECTORS  MB  AND  MA  ARE  ALSO  GENERATED. 

IMPLICIT  REAL*8 (A-H, 0- Z ) 

COMMON/ CM1A/ NLAT . NLONG 

COMMON/ CMS/ACZII.bTZI)  ,  v  ,  „  ,  %  ,  . 

COMMON  AG(ZB) ,Bg(ZC) ,GA(ZK) ,GB1(ZL) ,GB2(ZL) ,MA(ZM) , 

1MB  ( ZN  ) 

DIMENSION  BG( NLAT), AG (NLONG) ,GB1(2*NLAT-1) , 

GA / 3*NLONG -  3  J ,MA(N£0NG+1) ,MB(NLAT+2 ) 

LATX=NLAT+ 1 

FIND  BG  =  (ELEMENT  HEIGHT ) / 6 . 


g 

£ 

£ 

2 

2 

£ 

2 

FIND  AG 


f-l)+BG(J)) 


DO  10  J=1 .NLONG 
AG(Ji=A(jJ/3. 

GENERATE  GA,  ,  % 

GA( lj=2 . *(AG(1)+AG (NLONG ) ) 
DO  12  J= 2, NLONG 


DO  12  J=2,NL0 
K=2*LJ- 1 ) 
GaTk)=2*  (AG( 
GA(KUJ=AG(J- 
Kl=2*NLONG  , 
K2=3*NLONG-4 
DO  14  K=K1 ,K2 
GA(K)=0 . 


( 1) +AG( J) ) 


GA(K) =0 . 

GA( 3*NLONG - 3>= AG (NLONG ) 
GENERATE  DIRECTORY  VECTORS 
MB(1)=1 
D016  J=1.NLAT 

mb(j+u=2*j  ^ 

MB(NLAT+2)=2*(NLAT+1) 

MA ( 1 ) - 1 

DO  18  J=2. NLONG 
MA(JJ  =  2*(i-lL 

MACNLONG  ♦ 1 )  =  3*NLONG - 2 

RETURN 

END 
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Listing:  G0G3  EXEC 

ERASE  GAUSS2  *  Al 

COPY  GAUSS3  FORTRAN  Al  GAUSS2 

&STACK  CDIM 

& STACK  FILE 

X  GAUSS2  FORTRAN 

FORTGI  GAUSS2 

GLOBAL  TXTLIB  FORTMOD2  MOD2EEH 
FILEDEF  05  DISK 
LOAD  GAUSS2  (START 


Listing:  CDIM  XEDIT 

SET  CMS  TYPE  HT 
TOP 

*  ZB=NLONG 

c|zB/6/  *  * 

c|zC/3/  *  * 

*  Zl=NLONG*NLAT 
C|Z1/18 /  *  ** 

*  ZK=  3*NLONG- 3 

cJzK/iS/  *  * 

ZL=2*NLAT+1 

UZW7/ 

*  ZM=NLONG+l 

cJzM/7/  *  * 

*  ZN=NLAT+2 
C  |ZN/5/  *  * 

yi  syw®!(i^T"l) 

SET  CMSTYPE  RT 
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